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ABSTRACT 

The physical processes involved in the advective-acoustic instability are investigated with 2D nu- 
merical simulations. Simple toy models, developed in a companion paper, are used to describe the 
coupling between acoustic and entropy /vorticity waves, produced either by a stationary shock or by 
the deceleration of the flow. Using two Eulerian codes based on different second order upwind schemes, 
we confirm the results of the perturbative analysis. The numerical convergence with respect to the 
computation mesh size is studied with ID simulations. We demonstrate that the numerical accuracy 
of the quantities that depend on the physics of the shock is limited to a linear convergence. We 
argue that this property is likely to be true for most current numerical schemes dealing with SASI in 
the core-collapse problem, and could be solved by the use of advanced techniques for the numerical 
treatment of the shock. We propose a strategy to choose the mesh size for an accurate treatment of 
the advective-acoustic coupling in future numerical simulations. 

Subject headings: hydrodynamics — shock waves — - instabilities — - supernovae: general 



1. INTRODUCTION 

Most of our knowledge about the possible conse- 
quences of SASI on the core-collapse problem has 
been built, over the last 5 years, on the re- 
sults of multidimensiona l numer ical simulati ons (e.g. 
Blondin e t al. 2003; Scheck et all 2004; Burrows et al 
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; ^ondin & Mezzac appa I l2007i : iMarek fc Janka 
I; [fwak ami et al. 2008^^ Whether or not SASI 
contribute to overcome the explosion thresh- 
to kick the neutron star and alter its spin is 
debated. In addition to the fundamental un- 



certainties associated with the equation of state of 
dense matter or the numerical treatment of neu- 
trino transport, some difficulties a re simply related to 



multidimensional hy drodynamics (iBlondin et al 



2003; 



2007; 



phnishi et al. 2006: B londin fc Mezzacappa 112006 . 
[iwakami et al.i i2008i) . This latter difficulty is partly 
due to the complexity of the mechanism underlying 
SASI, which is at best unfamiliar, and possibly also 
affected by the different numerical techniques used by 
different groups. The present study aims at improv- 
ing our understanding of the instability mechanism at 
work by studying the advective-acoustic instability in 
the highly simplified set up introd uced in the first 
paper of this series (|Foglizzo 1 120081 . hereafter paper 
I). We note that a debate exi sts about the nature of 
this mechanism, as witnessed by Blondin & Mezzacappa 
(2006, h ereafter BM06), Fodizzo et al. (2007, hereafter 
FGS 07). iLamingl (120071) . lYamasaki &: Foglizzo I (|2008a ) 



and iLaming I ( 20081 ) . Thus we believe that a better un 



derstanding of the advective-acoustic instability in sim- 
ple examples can help recognise it in more complex sit- 
uations. The separation of the advective-acoustic cycle 
into two separate problems is necessary in order to iden- 
tify, between advected and acoustic perturbations, the 
conseque nces of each on e on the other, as s een on Fig. 7 
of iBlondin et all (2003") or Figs. 11-12 of IScheck et all 
(|2008l ). In paper I, the following questions were answered 
through a perturbative analysis: 



(i) what are the amplitudes of the entropy and vorticity 
waves generated by a shock perturbed by an acoustic 
wave propagating against the flow, towards the shock? 

(ii) what is the amplitude of the acoustic wave gen- 
erated by the deceleration of an entropy /vorticity wave 
through a localised gravitational potential? 

The first purpose of our study is thus to check the 
results of the perturbative analysis presented in paper I 
through numerical experiments, thus providing concrete 
examples of the coupling processes involved. 

The second purpose of this study is to gain confidence 
in the results of more elaborate numerical simulations 
by assessing their accuracy using our simple set up. The 
2D numerical simulations of BM06 showed some glob- 
ally good agreement with the perturbative analysis of 
FGSJ07. The typical error on the growth rate and the 
oscillation frequency of SASI, around 30%, was not small 
though. Could this be a concern for the many other sim- 
ulations which use a coarser mesh size? We wish to eval- 
uate quantitatively, using our simple toy model, to what 
extent the advective-acoustic instability can be affected 
by numerical resolution. 

The paper is organised as follows. In Sect. 2, the set up 
of the simulations is described and the numerical codes 
are presented. Sect. 3 illustrates qualitatively the two 
coupling processes involved in the advective-acoustic in- 
stability using 2D simulations. A quantitative analysis of 
these simulations is also performed which validates both 
the perturbative analysis and the numerical technique. 
In Sect. 4, we evaluate the rate of numerical convergence 
with respect to the mesh size, using a series of ID nu- 
merical simulations. While the accuracy of the acous- 
tic feedback produced by the flow gradients is quadratic 
with respect to the mesh size, the accuracy of the en- 
tropy wave produced by the shock depends on the mesh 
size only linearly. The linear phase of the full problem 
is simulated in Sect. 5, where the oscillation frequency 
and growth rate are compared to the results of the per- 
turbative analysis. The consequences of these numerical 
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Fig. 1. — Schematic view of the advective- acoustic cycle oc- 
curring in the toy model, separated in two sub problems. En- 
tropy /vorticity perturbations are noted as circular arrows, while 
acoustic waves are noted as wavy arrows. The linear coupling be- 
tween waves is measured by the efficiencies Qsh, Qv^ ^sh and 
7^v• 

difficulties for the simulations of core-collapse supernovae 
are discussed in Sect. 6. 

2. NUMERICAL TECHNIQUES AND SET UP OF THE 
SIMULATIONS 

2.1. Numerical techniques 

The gs:overni n^ equations a re solved using the AUS- 
MDV scheme (|Wada fc Liou hl994) , which is a second- 
order finite volume scheme. The former version of 
AUSMDV was called "advection ups tream splitting 
method" (AUSM) and developed by iLiou fc Stefienl 
(1993). AUSM is a remarkably simple upwind fiux vec- 
tor splitting scheme that treats the convective and pres- 
sure terms of the fiux function separately. In the AUS- 
MDV, a blending form of AUSM and fiux difference 
is used, and the robustness of AUSM in dealing with 
strong shocks is improved. A great advantage of this 
scheme is the reduction of numerical viscosity, which 
gives sharp preservation of fiuid interfaces and high res- 
olution fe ature as in the "piecew i se pa rabolic method" 
(PPM) of IColella fc Woodward I (|l984f ). Some advan- 
tages over PPM are simplicity and a lower computa- 
tional cost. In Sect. 4, the numerical results obtained 
with AUSMDV are compared with those computed us- 
ing RAMSES (|Tevssierll2002h . RAMSES is also a sec- 
ond order shock-capturing code. It uses the MUSCL- 
Hancok scheme to update the MHD equation. For 
the simulations presented in sect. 5, we used the Min- 
Mod slope limiter along; w ith the HLLD Riemann solver 
(|Mivoshi fc Kusano 1120051). which reduces to the HLLC 
Riemann solver (|Toro et al .1119941 ) in the hydrodynamic 
case dealt with in this paper. 

2.2. General set up 

In this section, we describe the problems we designed 
to illustrate the physical mechanisms underlying the 
advective-acoustic instability. Our "Problem 1" is aimed 
at studying the interaction of waves in a stationary sub- 
sonic fiow decelerated across a localised external po- 
tential, whereas "Problem 2" studies the interaction of 
waves with a stationary shock in a uniform potential. 
Both problems were described in detail in the linear 
approximation in paper I, and are schematically illus- 
trated by Fig. [TJ Let us recall that the stationary fiow 



is uniform in the x direction, and fiows along the z di- 
rection with a negative velocity. The ideal gas satisfies 
a polytropic equation of state with an adiabatic index 
7 = 4/3, and a measure of the entropy is defined as 
S = (log(p/p^))/(7 — 1). The horizontal size of the com- 
putation domain is noted L^. The index "1" refers to the 
supersonic fiow ahead of the shock {z > Zgh), and "in" 
refers to the subsonic region after the shock {z < Zsh)- 
A^in, vi and pi are determined by the Rankine-Hugoniot 
relations as follows: 
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/ 2 + (7-l)X^ 
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(1) 

(2) 
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where = — A^inCin- The incident Mach number is 
chosen as Mi = 5. Thus Min ^ 0.39. 

A region of deceleration extends over a width ~ 
centred on = 0, separating two uniform subsonic re- 
gions indexed by "in" and "out", respectively. The ex- 
ternal potential A^{z) responsible for the fiow gradients 
is defined by 
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The potential jump A<I> > is set by specifying the sound 
speed ratio cin/cout- 



A<l> = 



ML 
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Defining H = z^h — we adopt / H 



(5) 
0.1 and 



,2 



0.75 in this study, as in paper I. 



Time is normalised by Taac, which is a reference 
timescale associated to the advective-acoustic cycle de- 
fined as follow: 



I -Mi. 
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The advection time through the deceleration region ry 
is associated in paper I to a frequency cut-off c^v, above 
which the efficiency of acoustic feedback decreases: 



1 



TV 



(7) 
(8) 



Units are chosen such that cin = 1, Pin = 1 and H = 1. 
Since p = pc? and 7 = 4/3, then p\^ = 0.75 and Sin ~ 
—0.86. The reference timescale is thus r^ac ~ 4.2, and 
TV ~ 0.41, so that a;vTaac/27r ~ 1.6. 

Periodic boundary conditions are applied in the x- 
direction. Linear perturbations are characterised by 
their wavenumber kx = 27Tnx/Lx^ with = 4, and 
their frequency ujq. With this set of parameters, we 
expect from paper I a dominant mode Ux = 1 with a 
growth rate uJiTg^g^c =0.22 and an oscillation frequency 

UJrTaac/^TT = 1.13. 

With these parameters, the frequency c^ev below which 



3 



acoustic waves are evanescent in the z direction is 
^ev'^aac/27r = 0.96 in the uniform subsonic region be- 
fore deceleration, and a;g^*raac/27r = 1.20 after decelera- 
tion (Eq. (13) in paper I). For a;rTaac/27r = 1.13, acoustic 
waves are evanescent after the region of deceleration with 
an evanescence length ~ 1.9i^, deduced from Eq. (19) 
in paper I. 

2.3. Set up of 'Troblem r 

In "Problem 1", the flow is only composed of three 
parts, without a shock, and is thus entirely subsonic. 
Once the stationary unperturbed flow is well established 
on the computation grid, an entropy /vorticity wave is 
generated at the upper boundary, at z = 3. This wave 
is in pressure equilibrium {5p = 0). The corresponding 
perturbations of entropy SS and density Sp are defined 
as follows: 



SS = es cos {—coot + kxX + kzz) , 
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where = 10~^ is the parameter defining the amplitude 
of the entropy perturbation. The vertical wavenumber 
of an advected wave is kz = ujo/vin- The incompressible 
velocity perturbations Svx are Svz are chosen such that 
the vorticity Swy is the same as when produced by a 
shock (Eqs. (A6-A9) in paper I): 

k^cjocf^ SS 

^^^= ,,2 , — ' (11) 
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(12) 



(13) 



We choose free boundary conditions at the lower bound- 
ary {z = —5), sufficiently far from the shock to avoid any 
effect from a reflected wave. Between z = —2 and —5, we 
use an inhomogenous mesh whose interval increases grad- 
ually in the negative z-direction. We perform simulations 
with kx = 27t/Lx and different values of the frequency 
coq and mesh size Az. The results of the simulations are 
analysed in Sect. 3.1 and 3.2. 

2.4. Set up for 'Trohlem 2'' 

In our "Problem 2", the unperturbed stationary flow 
is composed of two semi-infinite uniform regions sepa- 
rated by a stationary shock. Once the steady flow is well 
established on the numerical grid, an acoustic wave is 
generated at the lower boundary of the computing box, 
at z = — 2 and propagates against the flow towards the 
shock. The density perturbation Sp^ the pressure per- 
turbation Sp and the velocity perturbations Svx and Svz 
are defined according to paper I as follows at the lower 
boundary: 
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X Cp COS (— c^o^ + kxX -\- k^ z) , (14) 

(15) 



Sv^ 



X Cp COS {—ujot + kxX -\- k^ z) ^ (16) 
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Fig. 2. — Production of an acoustic wave by the deceleration of 
a vorticity wave (Problem 1). The specific vorticity Swy/p (left) 
and the normalised pressure perturbation 6p/p (right) are shown at 
three successive times, before and after the advected wave reaches 
the deceleration region localised around z = (within the dashed 
lines). The par^rn^^rs are a;oTaac/27r = 2, and Ax = Az = 10""^. 

^ " COS (^—coot -\- kxX -\- k~ z) ^{17) 
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(18) 



Here = 10 ^ sets the amplitude of the density per- 
turbation. The vertical wavenumber k~ for an acoustic 
perturbation is given by Eq. (19) of paper I: 

Cin I - M 



2 • 



(19) 



We choose fixed boundary conditions at the upper 
boundary {z = 2). The results of the simulations are 
analysed in section 3.3 and 3.4. 

3. NUMERICAL ILLUSTRATION OF THE COUPLING 
PROCESSES AND COMPARISON WITH THE LINEAR 
ANALYSIS 

3.1. Acoustic feedback from the deceleration of a 
vorticity wave (Problem 1) 

The snapshots in Fig. [2] show the specific vorticity 
Swy/p (left column) and pressure perturbation Sp/p 
(right column) in the flow at three successive times, 
before and after the moment when the advected wave 
reaches the deceleration region. The right column of 
Fig. [2] demonstrates the absence of an acoustic perturba- 
tion until the advected wave reaches the region of deceler- 
ation. Two acoustic waves are then generated, propagat- 
ing upward and downward. This simple experiment gives 
a concrete illustration of the physical process described in 
analytical terms in paper I. In the bottom plots of Fig. [21 
the flow has reached the asymptotic regime described by 
a single frequency in paper I, in which a more quanti- 
tative comparison of coupling efficiencies can be made. 
Since the computation domain is finite, the numerical ex- 
periment is stopped before the acoustic waves reach the 
vertical boundaries of the computation box in order to 
avoid spurious reflections. The time needed to reach the 
asymptotic regime described by a single frequency in pa- 
per I depends strongly on the frequency of the wave, and 
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Fig. 3. — Efficiency (6pQ/pin)/^S of the production of acoustic 
waves by the deceleration of entropy /vorticity waves, measured at 
z = 0.5, as a function of ujq in Problem 1. The solid line shows 
the curve computed by a linear analysis (paper I). The results of 
numerical simulations are shown for different square mesh sizes 
Ax = Az = 5 X 10~^ (crosses), 2 x 10"^ (triangles) and 10"^ 



(circles). The results for Ax 
shown (pluses). 



2 X 10-2 Az = 10"2 are also 



can become prohibitively long close to the frequency of 
horizontal propagation cj^^. This can be understood by 
viewing the semi-infinite acoustic plane wave, involved in 
both Problems 1 and 2, as an infinite plane wave of fre- 
quency cjo, multiplied by a step function, whose Fourier 
transform involves a continuum of frequencies. In 1-D, 
all frequencies would propagate with the same velocity, 
and the shape of the wave packet would stay unchanged 
during propagation. In 2D however, the high frequency 
part of the acoustic spectrum co > coq propagates more 
vertically than the main component, while the low fre- 
quency part ujo > ^ > ^ev propagates more horizontally: 
this dispersion requires a longer numerical simulation, 
and thus a larger computational domain in order to avoid 
acoustic refiections. For this reason we have limited our 
investigation to the frequencies cJoTaac/^Tr = 1.5, 2, 4, 
and 6. Note that if the frequency of the perturbation had 
been chosen below the threshold of acoustic propagation 
{co < cjg^), the acoustic feedback would be evanescent 
above the deceleration region (paper I and Guilet, Sato 
& Foglizzo, in preparation). 

3.2. Measure of the acoustic feedback in Problem 1 

The amplitude of the acoustic feedback is measured in 
the numerical experiment by using a Fourier transform, 
in time, of the pressure perturbation over the period T = 
27r/uJo of the wave: 



(20) 



The symbols in Fig.[3]are measured at 2: = 0.5, in a region 
where the gravitational potential is uniform. The full line 
in Fig. [3] shows the expected efficiency {SpQ/pin)/SS of 
the acoustic feedback obtained by integrating the differ- 
ential system as in paper I. The good agreement with 
the perturbative calculation for a fine mesh (circles) con- 
firms the validity of both the perturbative calculation 
and the numerical code. Given the long horizontal wave- 
length of the perturbations, the results are insensitive to 
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Fig. 4. — Production of a vorticity wave by the interaction of 
an oblique acoustic wave with the shock (Problem 2). Swy/p (left) 
and 6p/p (right) are shown at three successive times, before and 
after the acoustic wave reaches the shock localised ai z = 1 (dashed 
line). A vorticity wave is generated and advected downward. The 
parameters are a;oTaac/27r = 2, and Ax = Az = 10~^ . 

an increase of the horizontal size Ax of the mesh (pluses 
and circles). As described in paper I, the efficiency of 
the acoustic feedback decreases for frequencies above the 
cut-off cjv ^ 1/tv- 

3.3. Entropy /vorticity produced by a shock perturbed by 
an acoustic wave (Problem 2) 

The upward propagation of the acoustic wave gener- 
ated at the lower boundary of the computation domain 
in Problem 2 is visible on the right column of Fig. [H 
The three snapshots illustrate the independence of ad- 
vected and acoustic perturbations in the uniform part of 
the ffow: the vorticity wave visible on the left column in 
Fig. [4] is generated only as the acoustic wave reaches the 
shock. This vorticity wave is then continuously generated 
by the shock and advected downward with the ffow. An 
entropy wave (not shown) is also generated at the shock, 
with the same appearance as the vorticity wave. The 
lower boundary condition in this experiment is chosen 
far enough so that the reffected acoustic wave generated 
at the shock does not have time to interact with the lower 
boundary. The efficiency of entropy /vorticity generation 
at the shock can be measured at the time correspond- 
ing to the bottom panel in Fig. [H and compared to the 
calculations of paper I. 

3.4. Measure of the entropy production in Problem 2 

According to Eqs. (30-31) of paper I, the amplitude 
5Sth. of the entropy wave produced by an acoustic wave 
reaching the shock is expected to be related to the fre- 
quency of the pressure wave as shown by the full line in 
Fig. El 
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Measuring the amplitude of the entropy wave produced 
by the shock in the numerical simulations is not straight- 
forward because of the presence of spurious high fre- 
quency oscillations, analysed in more details in the next 
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Fig. 5. — Dependence of SSo/(6p/pin), measured at z = 0.5, on 
the frequency ujq, in Problem 2. The sohd hue shows curve pre- 
dicted from hnear analysis (paper I). The result of numerical sim- 
ulations is shown for different mesh sizes Az = 2 x 10 ~^ (pluses), 
Az = 10~^ (squares), 5 x 10~^ (crosses), 2 x 10"*^ (triangles) and 
10~^ (circles) where Ax = 2 x 10~^. The filled points show the 
results for Ax = Az = 10~^. 

section. We choose to measure (at z = 0.5) its fun- 
damental Fourier component 5So at the frequency cjq, 
thus filtering out oscillations at higher frequency. The 
result is displayed in Fig. [5] for different frequencies and 
mesh sizes. We did not notice any dependence on the 
horizontal size Ax of the mesh, for the long horizontal 
wavelengths considered. The expectation of the pertur- 
bative calculation is confirmed, but the convergence to 
the analytical formula is apparently much slower than for 
Problem 1. The rate of convergence is analysed in the 
next section using ID simulations. 

4. ACCURACY OF THE NUMERICAL CONVERGENCE 

The dependence of the numerical error on the mesh size 
is easier to investigate using ID simulations because of 
the shorter computation time. Without excluding the 
possibility of additional difficulties in 2D, we demon- 
strate here that some numerical difficulties associated 
to the advective-acoustic coupling are already present in 
ID. The set up we use in this section is the same as used 
for the 2D simulations except that kx = 0. 

4.1. Quadratic convergence in Problem 1 

A series of numerical simulations of Problem 1 in ID 
with different mesh sizes and perturbation frequencies 
allowed us to measure the accuracy of the computation 
compared to the perturbative analysis as shown in Fig. [6] 
by the open squares. They are to be compared with 
the dotted line, whose slope of +2 illustrates second or- 
der convergence for this problem. Remembering that 
the accuracy of our numerical scheme is second order 
in space, it is satisfactory to find that the error dis- 
played in Fig. [6] is approximately quadratic with respect 
to the mesh size. The shortest wavelength in Problem 1 
is the wavelength 27ri'out/^o of advected perturbations 
after their deceleration, which is equal to ~ 0.12 for the 
frequency C(;oTaac/27r = 6. We conclude from Fig. [6] that 
our numerical treatment of advection, propagation and 
advective-acoustic coupling involved in Problem 1 is ac- 
curate at the percent level even when the shortest wave- 
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Fig. 6. — Numerical error as a function of the mesh size for 
Problem 1. The panels (a), (b) and (c) correspond to the cases of 
^o'?'aac/27r = 2, 4 and 6, respectively. The dotted lines, propor- 
tional to A^;^, illustrate the quadratic convergence. 

length is sampled by only ~ 10 grid zones. 

4.2. Linear convergence in Problem 2 

Applying the same test to Problem 2 is more compli- 
cated because of the high frequency oscillations already 
mentioned in Sect. 3. The shape of the entropy wave is 
shown in Fig. [71 for different frequencies and mesh sizes. 
The finer the mesh the higher the frequency of these spu- 
rious oscillations. We checked that the power involved in 
the Fourier component associated with these higher fre- 
quencies is always negligible compared to the main com- 
ponent. The Fourier component associated with the fun- 
damental frequency uoq converges slowly to the expected 
analytical value for a fine mesh. The squares in Fig. [8] 
show the numerical accuracy of the AUSMDV scheme 
for Problem 2, revealing a linear convergence with the 
mesh size (as shown by the dotted line of slope +1). We 
note that a coarse resolution can either underestimate 
or even overestimate the production of entropy at the 
shock. In order to show that this linear convergence is 
not a peculiarity of the AUSMDV scheme, these simula- 
tions were repeated with the code RAMSES. The results 
obtained with RAMSES are shown by the blacks circles 
in Fig. [8] (note that we also observed spurious high fre- 
quency oscillations in that case). They are comparable 
to those obtained using the AUSMDV scheme. Based 
on this comparison, we anticipate that all finite volume 
codes in which the treatment of the shock relies on an 
upwind technique are likely to share the same difficulty: 
quantities produced at the shock location, such as vor- 
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for the same three frequencies as in Fig. [6] The thick line, dotted 
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and thin lines correspond to the cases Az ■ 
respectively. 

ticity and entropy waves, or the reflected acoustic wave, 
are computed with a first order accuracy with respect to 
the mesh size. Likewise, we anticipate that all finite vol- 
ume codes will suffer from the presence of spurious high 
frequency oscillations similar to those described above. 
It is indeed well known that such codes are subject to 
this problem, especially in the case of standing shocks, 
as was reported by Colella & Woodward (1984). In the 
present case, the problem is made worse by the inter- 
action between the shock and the sound wave (in the 
absence of the latter, we barely detected high frequency 
oscillations, with an amplitude of the order of 0.5% of the 
am plitude of t he reflected entropy wave). As described 
by IColella fc Woodward \ (1984), any additional source 
of dissipation (artificial viscosity, grid translation) will 
result in a decrease of the amplitude of the oscillations. 
For example, with RA MSES, the u se of the Monotonised 
Central slope limiter (|Toro I [19971 ). which is known to 
be less dissipative than MinMod, resulted in the ampli- 
tude of the oscillations being about three times larger. 
However, the complete stabilisation of the oscillations 
(through the use of artificial viscosity for example) would 
most probably come at the cost of reducing the growth 
rate, which we show in Sect. 5 not to be affected by the 
oscillations. 

5. EIGENFREQUENCY IN THE FULL TOY MODEL 

The full toy model has been simulated in order to mea- 
sure the oscillation frequency uJr and growth rate uJi of 
the dominant eigenmode for A^i = 5, H\/ / H = 0.1, 
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Fig. 8. — Numerical error of the quantity \6So/6Sth\ as a func- 
tion of the mesh size for Problem 2. The frequencies are the same as 
in Fig. [6] The empty squares and filled circles were obtained with 
the AUSMDV scheme and the code RAMSES respectively. The 
dotted lines, proportional to Az, illustrate the linear convergence. 



simulation is the numerical relaxation of the unperturbed 
flow on the computational grid, which can result in a slow 
drift of the shock. The stationary flow is constructed by 
first obtaining a stationary subsonic flow in the gravi- 
tational potential, and then choose the upstream flow 
such that a shock is stationary at 2: = 2:sh- As a result 
of numerical discretization, the upstream mach number 
may slightly differ from A^i = 5, by a few per cents. This 
difference is taken into account in the perturbative calcu- 
lation of the reference eigenfrequency. Perturbations are 
incorporated as a random noise in the transverse veloc- 
ity at the level of 10% of the flow velocity in the uniform 
region between z = 0.3 and z — 0.9. The linear evolution 
is dominated by the mode n^c = 1, as expected from the 
linear stability analysis. The comparison with the per- 
turbative calculation is shown in Fig. [9l The oscillation 
frequency and growth rate, determined numerically, are 
accurate to about 5% for Az < 10~^, suggesting that the 
spurious high frequency oscillations revealed in Sect. 3.4 
and 4.2 have a minor effect on the eigenfrequency of the 
most unstable mode. The slight excess of the growth rate 
^i,sim in Fig. [9] may be related to the fact that entropy 
and vorticity pertrubations are slightly overproduced at 
the shock, as seen in Fig. 5 for Problem 2. This effect, 
however, should be partially compensated by the slight 
underproduction of the acoustic feedback in Problem 1 
(Fig. 3). 

A significant damping of the instability (~ 14%) occurs if 
the grid is too coarse (A2: = 2 x 10~^) but even then, the 
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Fig. 9. — Growth rate uji^sim and oscillation frequency a;r,sim 
of the most unstable mode (ux = 1) measured in a numerical 
simulation of the full toy model, compared to the values o^^^th? ^r,th 
obtained from the perturbative analysis (paper I). The parameters 
are L^/H = 4, H^/H = 0.1, Mi = 5, cfjc^^^ = 0.75, and 
Ax = 10~^. Perturbations were initiated with a random noise. 
Error bars are associated to the fitting procedure. 

oscillation frequency is accurate within 5%. The surpris- 
ing accuracy of the oscillation frequency can be under- 
stood by the fact that the oscillation timescale is closely 
related to the timescale, for an advective- acoustic cycle 
between the shock and the deceleration region. Since the 
position of the acoustic feedback is set by the external 
potential in our toy model, this timescale barely depends 
on the numerical resolution. One must keep in mind that 
in a realistic flow where gradients are due to cooling pro- 
cesses, a change of numerical resolution could influence 
the position of the deceleration region, and could thus 
affect the oscillation timescale of the instability. 

6. CONSEQUENCES FOR CORE-COLLAPSE SIMULATIONS 

The results of our numerical experiments can be help- 
ful to choose the mesh size in future simulations of a 
collapsing stellar core, both at the shock and near the 
neutron star, in order to make sure that the physics of 
SASI is correctly treated, at least in the linear regime. 
Of course, the influence of SASI on the mechanism of 
core-collapse supernovae depends on non-linear quanti- 
ties such as the amplitude of the shock oscillations, the 
advection time through the gain region, or the spectral 
distribution of energy below the shock. Characterising 
which of the non-linear properties of SASI are most sen- 
sitive to the numerical technique is beyond the scope of 
the present study, and will be investigated in a forthcom- 
ing publication. We believe however that the coupling 
between entropy, vorticity and pressure is likely to play 
an important role even in the non linear regime of SASI, 



both through the flow gradients and at the shock. The 
wide range of frequencies involved in the non linear evolu- 
tion of SASI (e.g. Yoshida et al. 2007) suggests that the 
accuracy of the numerical treatment should not be lim- 
ited to the low frequency of the most unstable mode. In 
this sense, the numerical constraints deduced from our 
linear analysis should be considered as a minimum re- 
quirement, even-though some non-linear consequences of 
SASI may be less sensitive to numerical resolution than 
others: the addition of numerical errors with opposite 
signs, mentioned in Sect. 5, may contribute to the com- 
plex, non monotonic dependence of the explosion time 
with respect to the numerical resolution, observed by 
Murphy & Burrows (2008). 

6.1. Mesh size in the deceleration region 

When the shock stalls above the proto-neutron star, 
the flow deceleration close to the neutron star is domi- 
nated by cooling processes much more than by gravity, 
and the advective-acoustic coupling there is not adia- 
batic. By making the choice of simplicity, our toy model 
does not aim at reproducing quantitatively the efficiency 
of the acoustic feedback in a non-adiabatic flow. It helps 
understand that a simulation with a coarse grid in the 
vicinity of the neutron star may be unable to take into 
account a possible acoustic feedback from this region, 
simply because advected perturbations are numerically 
damped before reaching it. Let us consider a numeri- 
cal simulation of an advective-acoustic cycle dominated 
by the oscillation frequency ujq- The choice of the mesh 
size close to the surface of the neutron star is not ob- 
vious because the wavelength of advected perturbations 
Aadv ^ 27rv{r)/uJo shrinks as the gas is decelerated. For- 
tunately, an accurate advection of this perturbation is 
needed only down to the region where most of the acous- 
tic feedback is generated, adiabatic or not. Since the 
timescale of the advective-acoustic cycle is larger than 
the advection timescale, and comparable to the oscilla- 
tion timescale 2i\ juji of the fundamental mode, the region 
of feedback is necessarily above the radius rin reached by 
the gas during one SASI oscillation. According to Figs. 4 
and 5 of FGS JOT, the dominant mode is the fundamental 
one {ujf = uo) if the shock is close to the neutron star, 
or the first harmonic {oOf ^ 2cc;o) if the shock distance is 
large enough, rin is thus defined by: 



rsh ^ 27r 

1^1 ^f* 



(22) 



A possible strategy to choose the mesh size Ann in the 
inner region of the flow could be to make sure that the 
advected perturbations are correctly advected down to 
this radius rin- Denoting by N the number of grid zones 
per wavelength required for an accurate advection and 
acoustic coupling of vorticity perturbations, the maximal 
mesh size Arin near the radius rin should be 



Arin = ^ — v(rin . 



(23) 



Our illustration in Fig. [6] suggests N ~ 10. Of course, the 
precise value of N depends on the numerical technique 
used and is expected to vary from code to code but is 
likely to remain of the same order as our estimate. In 
any case, Eq. ([23|) will be useful for future numerical 
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Fig. 10. — Two-dimensional distribution of the power in the a 
fundamental mode, |5S'o/5S'th|^, function of A(^/Az and X^dv/^^ 
obtained in the ID simulations of problem 2. 

simulations involving SASI, as a consistency check that 
the advective-acoustic feedback is properly resolved, at 
least for the fundamental mode. 

6.2. Mesh size near the stalled shock 

Our study of Problem 2 has identified the difficulty 
of accurately calculating the entropy generated by the 
shock in a numerical simulations. This difficulty is likely 
to affect any physical quantity depending on the physics 
of the shock, such as the vorticity and the amplitude of 
reffected pressure waves. In this sense, all the numeri- 
cal simulations of core-collapse involving SASI must face 
a similar difficulty with the numerical treatment of the 
shock. 

We argue that this difficulty is not speciffc to the lin- 
ear regime of the instability. In the non linear regime 
of SASI, as long as the shock continues to play a fun- 
damental role by generating entropy and vorticity per- 
turbations, the accuracy of the quantities depending on 
its behaviour are likely to be affected by this first or- 
der convergence. However, the details and precise conse- 
quences of this issue in that case remain an open issue at 
the present time. Answering these questions will require 
more realistic simulations, coupling both problems and 
carried to the non linear regime. 

Should the grid size be able to resolve the displace- 
ment of the shock for a better accuracy ? According to 
the perturbative analysis, the shock displacement A( is 
related to the entropy perturbation 5S by Eq. (16) of 
paper I: 



AC: 



ss 



1 



(24) 



We show on Fig. [TO] the accuracy of the numerical sim- 
ulation, compared to the linear calculation, depending 
on how the grid sizes compares to both the shock dis- 
placement A( and the advection wavelength Aadv = 
27r I Vin I /cJo 7 in ID calculations. Non linear effects be- 
come dominant for > Aadv /1 00- In the linear regime 
(AC < Aadv/100), an accuracy of 10% requires Az < 
Aadv/100. Resolving the shock displacement does not 
seem to be a crucial condition for the computation of 
the entropy production. 

Since the exact properties of numerical convergence 
vary from a numerical scheme to another, it is not possi- 



ble here to determine the real accuracy of existing numer- 
ical simulations involving SASI. At best we can estimate 
what would be the accuracy of our AUSMDV scheme in 
the conditions used by various authors. The mesh size 
Arsh at the radius of the stalled shock in published sim- 
ulations varies depending on their complexity and the 
size of their outer boundary. We es timated Arsh 1 km 
in the 2D simul ations of BMQ6 and l Schec k et al? ("20 081) . 
Arsh 2 km in lOhnishi et al.| (l2006h and Iwakami et al.l 
(2008'), and Argh - 5 km in Burro ws et al.t (2006). Es- 
timating the value of the ratio Aadv/Argh is possible by 
identifying ujq with the oscillation frequency of the domi- 
nant mode. We esti mated Aadv/Argh ~ 200 in BM06 and 
IScheck et all (j2008[ ). which seems marginally sufficient to 
obtain a 10% accuracy from the point of view of Fig. [TOl 
The discrepancy of 30%, noted by FGSJ07 between the 
numerical results of BM06 and the perturbative analysis 
when the shock distance increases, may be related to the 
fact that the instability becomes dominated by the first 
harmonic rather than the fundamental mode. The corre- 
spondingly deeper coupling region may require a smaller 
mesh size, as already noted in FGSJ07 on the basis of 
the structure of the eigenfunction. Remembering that 
the mesh size in BM06 is one of the finest among the ex- 
isting core-collapse simulations, particular attention on 
this issue seems necessary for the future simulations in 
which SASI could play an important role. 

7. CONCLUSIONS 

• A toy model has been used to illustrate through 
numerical experiments the coupling processes de- 
scribed in mathematical terms in paper I. Despite 
the high degree of simplification of our toy model, 
in particular the adiabatic hypothesis and the very 
local character of the deceleration region, these 
simulations can help us build our intuition about 
the physics of the advective-acoustic instability and 
better recognise it when present in numerical sim- 
ulations. 

• The results of the perturbative approach have been 
confirmed quantitatively by our numerical simula- 
tions. 

• We have studied the effect of the mesh size on the 
accuracy of the numerical calculation. This will 
prove useful in the future to improve the reliability 
of the hydrodynamical part of simulations involv- 
ing SASI in the core-collapse problem. We have 
proposed a conservative estimate of the desired 
mesh size close to the neutron star, which guar- 
antees that the dominant acoustic feedback from 
advected perturbations is correctly taken into ac- 
count. 

• The difficulties associated with the numerical treat- 
ment of the shock have direct consequences on the 
accuracy with which the ffow resulting from SASI is 
calculated: without a special numerical effort, the 
convergence of the computation of the growth time 
and oscillation frequency of SASI is reduced to first 
order even if the numerical scheme converges with a 
higher order away from the shock. Among the pub- 
lished simulations of SASI, only the 2D simulations 
with the finest grid seem to be able to estimate the 
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entropy and vorticity production at the shock with 
a < 10% accuracy. The importance of an accurate 
treatment of S ASI in the core-cohapse problem may 
make it worth implementing advanced techniques 
for the numerical treatment of the shock in future 
sim ulations, such as the level s et method for exam- 
ple (jSethian fc Smereka |[2QQ3[ ). 
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